BST213: Final exam
Install needed packages
# install.packages("gdata")
library(gdata)
## gdata: Unable to locate valid perl interpreter
## gdata:
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata:
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
##
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
##
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
##
## Attaching package: 'gdata'
## The following object is masked from 'package:stats':
##
## nobs
## The following object is masked from 'package:utils':
##
## object.size
## The following object is masked from 'package:base':
##
## startsWith
# install.packages("car")
library(car)
## Warning: package 'car' was built under R version 3.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.5.2
# install.packages("lubridate")
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
# install.packages("gmodels")
library(gmodels)
# install.packages("RQuantLib")
library(RQuantLib)
## Warning: package 'RQuantLib' was built under R version 3.5.3
# install.packages("e1071")
library(e1071)
## Warning: package 'e1071' was built under R version 3.5.3
# install.packages("ggplot2")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
# intall.packages("plotly")
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.3
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# install.packages("survival")
library(survival)
Load database
# Load database
pSERG <- read.csv("D:\\BST213 APPLIED REGRESSION FOR CLINICAL RESEARCH\\FINAL EXAM\\pSERG.csv")
Data cleaning and data transformation
# Transform date of status epilepticus into date format
pSERG$DATESEIZURE <- as.Date(pSERG$DATESEIZURE, format = "%m/%d/%Y")
# Order by date of status epilepticus
pSERG <- pSERG[order(pSERG$PATIENT_LABEL, pSERG$DATESEIZURE), ]
# Delete duplicate episodes from the same patient
pSERG <- pSERG[!duplicated(pSERG$PATIENT_LABEL), ]
# Delete patients with unknown age
pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS <- as.numeric(as.character(pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS))
## Warning: NAs introduced by coercion
pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS <- as.numeric(as.character(pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS))
## Warning: NAs introduced by coercion
pSERG$ageyears <- pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS + (pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS/12)
pSERG <- pSERG[complete.cases(pSERG[ ,"ageyears"]), ]
# Delete patients with unknown sex
pSERG <- pSERG[which(pSERG$SEX == "male" | pSERG$SEX == "female"), ]
pSERG$SEX <- droplevels(pSERG$SEX)
# Delete patients with control (non-refractory) SE
pSERG <- pSERG[which(pSERG$SE_GROUP == "refractory_case"), ]
pSERG$SE_GROUP <- droplevels(pSERG$SE_GROUP)
# Delete patients with unknown hospital onset
pSERG <- pSERG[which(pSERG$HOSPITALONSET == "no" | pSERG$HOSPITALONSET == "yes"), ]
pSERG$HOSPITALONSET <- droplevels(pSERG$HOSPITALONSET)
# Transform BZDTIME.0 to numeric
pSERG$BZDTIME.0 <- as.numeric(as.character(pSERG$BZDTIME.0))
## Warning: NAs introduced by coercion
# Delete patients with unknown time to first BZD
pSERG <- pSERG[complete.cases(pSERG[ , "BZDTIME.0"]), ]
# Transform AEDTIME.0 to numeric
pSERG$AEDTIME.0 <- as.numeric(as.character(pSERG$AEDTIME.0))
## Warning: NAs introduced by coercion
# Delete patients with unknown time to first non-BZD-AED
pSERG <- pSERG[complete.cases(pSERG[ , "AEDTIME.0"]), ]
# Delete patients with unknown type of SE (continuous vs intermittent)
pSERG <- pSERG[which(pSERG$TYPESTATUS == "continuous" | pSERG$TYPESTATUS == "intermittent"), ]
pSERG$TYPESTATUS <- droplevels(pSERG$TYPESTATUS)
# Create convulsive duration in minutes and eliminate patients with unknown convulsive duration
pSERG$CONVULSIVEDURATION <- as.numeric(as.character(pSERG$CONVULSIVEDURATION))
## Warning: NAs introduced by coercion
pSERG$CONVULSIVEmin <- pSERG$CONVULSIVEDURATION
pSERG$CONVULSIVEhr <- pSERG$CONVULSIVEDURATION * 60
pSERG$convulsivedurationmin <- ifelse(pSERG$CONVULSIVEDURATIONUNITS=="min", pSERG$CONVULSIVEmin, pSERG$CONVULSIVEhr)
pSERG <- pSERG[complete.cases(pSERG[ , "convulsivedurationmin"]), ]
# Delete patients with unknown race
pSERG <- pSERG[complete.cases(pSERG[ , "RACE"]), ]
pSERG$RACE <- droplevels(pSERG$RACE)
# Delete patients with unknown or nonsensical time of SE
pSERG$TIMESEIZURE_HOURS <- as.numeric(as.character(pSERG$TIMESEIZURE_HOURS))
## Warning: NAs introduced by coercion
pSERG <- pSERG[complete.cases(pSERG[ ,"TIMESEIZURE_HOURS"]), ]
pSERG$TIMESEIZURE_HOURS <- ifelse(pSERG$TIMESEIZURE_HOURS >= 100, pSERG$TIMESEIZURE_HOURS/100, pSERG$TIMESEIZURE_HOURS)
# Create variable day/night
pSERG$day <- ifelse(pSERG$TIMESEIZURE_HOURS >= 8 & pSERG$TIMESEIZURE_HOURS < 20, 1, 0)
# Delete patients with time of SE onset
pSERG <- pSERG[complete.cases(pSERG[ , "day"]), ]
###############VARIABLE CREATION#########################
# Divide race into White and non-white
pSERG$white <- ifelse(pSERG$RACE == "white", 1, 0)
# Create variable delay
pSERG$delay[grepl("delay", pSERG$PAST)] <- 1
pSERG$delay[!grepl("delay", pSERG$PAST)] <- 0
# Create variable cerebral palsy
pSERG$palsy[grepl("palsy", pSERG$PAST)] <- 1
pSERG$palsy[!grepl("palsy", pSERG$PAST)] <- 0
# Create variable cerebral palsy or delay
pSERG$delay_or_cerebralpalsy <- ifelse(pSERG$delay==1 | pSERG$palsy==1, 1, 0)
# Create variable none (no neurological comorbidities)
pSERG$none[grepl("none", pSERG$PAST)] <- 1
pSERG$none[!grepl("none", pSERG$PAST)] <- 0
# Create variable prior epilepsy
pSERG$priorepilepsy[grepl("epi",pSERG$PAST)] <- 1
pSERG$priorepilepsy[!grepl("epi",pSERG$PAST)] <- 0
# Create variable prior status
pSERG$status[grepl("status",pSERG$PAST)] <- 1
pSERG$status[!grepl("status",pSERG$PAST)] <- 0
# Create variable of at least one continuous infusion
pSERG$CI <- ifelse(!(pSERG$CONTMED.0==""), 1, 0)
# Transform CI time into numeric
pSERG$CONTTIME.0 <- as.numeric(as.character(pSERG$CONTTIME.0))
## Warning: NAs introduced by coercion
# Create ICU stay in days
pSERG$ICU_DURATION <- as.numeric(as.character(pSERG$ICU_DURATION))
## Warning: NAs introduced by coercion
pSERG$ICUhours <- pSERG$ICU_DURATION/24
pSERG$ICUdays <- pSERG$ICU_DURATION
pSERG$ICUdurationdays <- ifelse(pSERG$ICU_UNITS=="days", pSERG$ICUdays, pSERG$ICUhours)
# Transform EMS arrival to numeric
pSERG$EMSARRIVAL <- as.numeric(as.character(pSERG$EMSARRIVAL))
## Warning: NAs introduced by coercion
# Transform time to hospital arrival to numeric
pSERG$HOSPITALARRIVAL <- as.numeric(as.character(pSERG$HOSPITALARRIVAL))
## Warning: NAs introduced by coercion
# Create variable BZD before hospital arrival
pSERG$BZDbeforehospital <- ifelse(pSERG$BZDTIME.0 < pSERG$HOSPITALARRIVAL, 1, 0)
# Reclasify etiology
pSERG$etiology2 <- recode(pSERG$ETIOLOGY,
"'genetic' = 'genetic';
'metabolic'= 'metabolic';
'other' = 'other';
'structural' = 'structural';
'unknown' = 'unknown';
'' = 'unknown'")
# Structural etiology
pSERG$structuraletiology <- ifelse(pSERG$etiology2 == "structural", 1, 0)
# Create variable early in the academic year
pSERG$dateSE <- as.POSIXct(pSERG$DATESEIZURE, format = "%m/%d/%Y")
pSERG$monthSE <- month(pSERG$dateSE, label = FALSE)
pSERG$earlyacademicyear <- ifelse(pSERG$monthSE >= 7 & pSERG$monthSE <= 12, 1, 0)
# Create variable awareness of delays in time to treatment
pSERG$yearSE <- year(pSERG$dateSE)
pSERG$awareness <- ifelse(pSERG$yearSE >= 2015, 1, 0)
# Number of days between the beginning and end of the study
days_of_study <- seq(as.Date("2011-6-1"), as.Date("2019-5-1"), by="days")
# Proportion of holidays in all the days of the study
holiday_days_of_study <- as.factor(ifelse(isBusinessDay(calendar = "UnitedStates", dates = days_of_study), 0, 1))
# Create variable weekends and holidays
pSERG$holiday <- as.factor(ifelse(isBusinessDay(calendar = "UnitedStates", dates = pSERG$DATESEIZURE), 0, 1))
Demographics
# Proportion of patients with ID
CrossTable(pSERG$delay_or_cerebralpalsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 149 | 151 |
## | 0.497 | 0.503 |
## |-----------|-----------|
##
##
##
##
# Age
nobs(pSERG$ageyears)
## [1] 300
summary(pSERG$ageyears)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08333 1.25000 3.83333 5.70721 9.25000 20.74167
sd(pSERG$ageyears)
## [1] 5.171068
# Age in patients with ID
nobs(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
## [1] 151
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08333 2.20833 4.91667 6.37369 9.50000 20.74167
sd(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
## [1] 5.069913
# Age in patients without ID
nobs(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
## [1] 149
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.125 1.000 2.500 5.032 8.917 19.306
sd(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
## [1] 5.201736
# Sex
CrossTable(pSERG$SEX)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | female | male |
## |-----------|-----------|
## | 128 | 172 |
## | 0.427 | 0.573 |
## |-----------|-----------|
##
##
##
##
# Sex in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$SEX)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | female | male |
## |-----------|-----------|
## | 68 | 83 |
## | 0.450 | 0.550 |
## |-----------|-----------|
##
##
##
##
# Sex in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$SEX)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | female | male |
## |-----------|-----------|
## | 60 | 89 |
## | 0.403 | 0.597 |
## |-----------|-----------|
##
##
##
##
# Race
CrossTable(pSERG$white)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 110 | 190 |
## | 0.367 | 0.633 |
## |-----------|-----------|
##
##
##
##
# Race in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$white)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 52 | 99 |
## | 0.344 | 0.656 |
## |-----------|-----------|
##
##
##
##
# Race in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$white)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 58 | 91 |
## | 0.389 | 0.611 |
## |-----------|-----------|
##
##
##
##
# Prior epilepsy
CrossTable(pSERG$priorepilepsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 156 | 144 |
## | 0.520 | 0.480 |
## |-----------|-----------|
##
##
##
##
# Prior epilepsy in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$priorepilepsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 40 | 111 |
## | 0.265 | 0.735 |
## |-----------|-----------|
##
##
##
##
# Prior epilepsy in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$priorepilepsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 116 | 33 |
## | 0.779 | 0.221 |
## |-----------|-----------|
##
##
##
##
# Structural etiology
CrossTable(pSERG$structuraletiology)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 222 | 78 |
## | 0.740 | 0.260 |
## |-----------|-----------|
##
##
##
##
# Etiology in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$structuraletiology)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 111 | 40 |
## | 0.735 | 0.265 |
## |-----------|-----------|
##
##
##
##
# Etiology in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$structuraletiology)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 111 | 38 |
## | 0.745 | 0.255 |
## |-----------|-----------|
##
##
##
##
# Hospital onset
CrossTable(pSERG$HOSPITALONSET)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | no | yes |
## |-----------|-----------|
## | 207 | 93 |
## | 0.690 | 0.310 |
## |-----------|-----------|
##
##
##
##
# Hospital onset in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$HOSPITALONSET)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | no | yes |
## |-----------|-----------|
## | 107 | 44 |
## | 0.709 | 0.291 |
## |-----------|-----------|
##
##
##
##
# Hospital onset in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$HOSPITALONSET)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | no | yes |
## |-----------|-----------|
## | 100 | 49 |
## | 0.671 | 0.329 |
## |-----------|-----------|
##
##
##
##
# Type of status epilepticus
CrossTable(pSERG$TYPESTATUS)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | continuous | intermittent |
## |--------------|--------------|
## | 97 | 203 |
## | 0.323 | 0.677 |
## |--------------|--------------|
##
##
##
##
# Type of SE in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$TYPESTATUS)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | continuous | intermittent |
## |--------------|--------------|
## | 52 | 99 |
## | 0.344 | 0.656 |
## |--------------|--------------|
##
##
##
##
# Type of SE in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$TYPESTATUS)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | continuous | intermittent |
## |--------------|--------------|
## | 45 | 104 |
## | 0.302 | 0.698 |
## |--------------|--------------|
##
##
##
##
# Day versus night
CrossTable(pSERG$day)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 126 | 174 |
## | 0.420 | 0.580 |
## |-----------|-----------|
##
##
##
##
# Day versus night in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$day)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 71 | 80 |
## | 0.470 | 0.530 |
## |-----------|-----------|
##
##
##
##
# Day versus night in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$day)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 55 | 94 |
## | 0.369 | 0.631 |
## |-----------|-----------|
##
##
##
##
# Weekday versus weekend/holiday
CrossTable(pSERG$holiday)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 208 | 92 |
## | 0.693 | 0.307 |
## |-----------|-----------|
##
##
##
##
# Weekday versus weekend/holiday in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$holiday)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 109 | 42 |
## | 0.722 | 0.278 |
## |-----------|-----------|
##
##
##
##
# Weekday versus weekend/holiday in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$holiday)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 99 | 50 |
## | 0.664 | 0.336 |
## |-----------|-----------|
##
##
##
##
# Season
CrossTable(pSERG$earlyacademicyear)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 160 | 140 |
## | 0.533 | 0.467 |
## |-----------|-----------|
##
##
##
##
# Season in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$earlyacademicyear)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 77 | 74 |
## | 0.510 | 0.490 |
## |-----------|-----------|
##
##
##
##
# Season in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$earlyacademicyear)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 83 | 66 |
## | 0.557 | 0.443 |
## |-----------|-----------|
##
##
##
##
# Convulsive duration
nobs(pSERG$convulsivedurationmin)
## [1] 300
summary(pSERG$convulsivedurationmin)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 60.0 124.5 2184.0 281.5 172800.0
# Convulsive duration in patients with ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$convulsivedurationmin)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 75.0 140.0 984.6 300.0 37440.0
# Convulsive duration in patients without ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$convulsivedurationmin)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 53 120 3400 245 172800
# Length of stay
nobs(pSERG$ICUdurationdays)
## [1] 287
summary(pSERG$ICUdurationdays)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 2.00 4.00 12.08 10.94 180.00 13
# Length of stay in patients with ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ICUdurationdays)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 2.000 4.000 8.532 9.750 73.000 6
# Length of stay in patients without ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ICUdurationdays)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 2.000 4.083 15.695 12.000 180.000 7
# Mortality
CrossTable(pSERG$ALIVE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 300
##
##
## | | no | yes |
## |-----------|-----------|-----------|
## | 2 | 10 | 288 |
## | 0.007 | 0.033 | 0.960 |
## |-----------|-----------|-----------|
##
##
##
##
# Febrile in patients in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ALIVE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 151
##
##
## | | no | yes |
## |-----------|-----------|-----------|
## | 1 | 3 | 147 |
## | 0.007 | 0.020 | 0.974 |
## |-----------|-----------|-----------|
##
##
##
##
# Febrile in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ALIVE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 149
##
##
## | | no | yes |
## |-----------|-----------|-----------|
## | 1 | 7 | 141 |
## | 0.007 | 0.047 | 0.946 |
## |-----------|-----------|-----------|
##
##
##
##
Univariate analysis time to first BZD
## Outcome: time to first treatment
summary(pSERG$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 17.00 65.53 45.00 1440.00
sd(pSERG$BZDTIME.0)
## [1] 162.3031
skewness(pSERG$BZDTIME.0)
## [1] 5.255733
kurtosis(pSERG$BZDTIME.0)
## [1] 33.28925
plot_firstBZD <- ggplot(aes(x=BZDTIME.0), data=pSERG) + geom_density(color="green", fill="green") + geom_vline(xintercept=mean(pSERG$BZDTIME.0), color="red", lwd=1) + geom_vline(xintercept=median(pSERG$BZDTIME.0), color="blue", lwd=1) + theme_classic()
plot_firstBZD
ggplotly(plot_firstBZD)
# Time to first treatment
survfit(Surv(pSERG$BZDTIME.0) ~ 1)
## Call: survfit(formula = Surv(pSERG$BZDTIME.0) ~ 1)
##
## n events median 0.95LCL 0.95UCL
## 300 300 17 14 20
plot(survfit(Surv(pSERG$BZDTIME.0) ~ 1), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("purple4"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
# Time to first BZD depending on ID
summary(pSERG[which(pSERG$delay_or_cerebralpalsy == 1), ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.00 61.91 48.50 1264.00
summary(pSERG[which(pSERG$delay_or_cerebralpalsy == 0), ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 5.0 20.0 69.2 42.0 1440.0
survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)
## Call:
## survdiff(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## pSERG$delay_or_cerebralpalsy=0 149 76.8 78.1 0.0208 0.0649
## pSERG$delay_or_cerebralpalsy=1 151 79.7 78.4 0.0207 0.0649
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.8
pchisq(survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)$chisq, length(survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)$n)-1, lower.tail = FALSE)
## [1] 0.7988799
# Figure time to first BZD by ID
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("ID", "No ID"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by age
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$ageyears))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$ageyears)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$ageyears 0.0004671 1.0004672 0.0108028 0.043 0.966
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$ageyears 1 0.9995 0.9795 1.022
##
## Concordance= 0.518 (se = 0.021 )
## Rsquare= 0 (max possible= 1 )
## Likelihood ratio test= 0 on 1 df, p=1
## Wald test = 0 on 1 df, p=1
## Score (logrank) test = 0 on 1 df, p=1
# Time to first treatment by sex
summary(pSERG[pSERG$SEX=="male", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 20.00 54.38 45.00 720.00
sd(pSERG[pSERG$SEX=="male", ]$BZDTIME.0)
## [1] 107.4961
summary(pSERG[pSERG$SEX=="female", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.00 80.52 46.00 1440.00
sd(pSERG[pSERG$SEX=="female", ]$BZDTIME.0)
## [1] 214.6275
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$SEX))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$SEX)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$SEXmale 0.04609 1.04717 0.11812 0.39 0.696
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$SEXmale 1.047 0.955 0.8308 1.32
##
## Concordance= 0.491 (se = 0.018 )
## Rsquare= 0.001 (max possible= 1 )
## Likelihood ratio test= 0.15 on 1 df, p=0.7
## Wald test = 0.15 on 1 df, p=0.7
## Score (logrank) test = 0.15 on 1 df, p=0.7
# Figure time to first BZD by sex
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$SEX), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=levels(pSERG$SEX),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by race
summary(pSERG[pSERG$white==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 19.00 67.62 47.75 1440.00
sd(pSERG[pSERG$white==1, ]$BZDTIME.0)
## [1] 157.3164
summary(pSERG[pSERG$white==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.00 61.93 39.50 1264.00
sd(pSERG[pSERG$white==0, ]$BZDTIME.0)
## [1] 171.2509
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$white))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$white)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$white -0.1113 0.8947 0.1203 -0.925 0.355
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$white 0.8947 1.118 0.7068 1.133
##
## Concordance= 0.516 (se = 0.017 )
## Rsquare= 0.003 (max possible= 1 )
## Likelihood ratio test= 0.85 on 1 df, p=0.4
## Wald test = 0.86 on 1 df, p=0.4
## Score (logrank) test = 0.86 on 1 df, p=0.4
# Figure time to first BZD by race
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$white), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("non-white", "white"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by prior epilepsy
summary(pSERG[pSERG$priorepilepsy==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.00 59.78 51.25 660.00
sd(pSERG[pSERG$priorepilepsy==1, ]$BZDTIME.0)
## [1] 119.3409
summary(pSERG[pSERG$priorepilepsy==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.00 20.00 70.85 40.50 1440.00
sd(pSERG[pSERG$priorepilepsy==0, ]$BZDTIME.0)
## [1] 193.9492
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$priorepilepsy 0.0274 1.0278 0.1166 0.235 0.814
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$priorepilepsy 1.028 0.973 0.8177 1.292
##
## Concordance= 0.511 (se = 0.018 )
## Rsquare= 0 (max possible= 1 )
## Likelihood ratio test= 0.06 on 1 df, p=0.8
## Wald test = 0.06 on 1 df, p=0.8
## Score (logrank) test = 0.06 on 1 df, p=0.8
# Figure time to first BZD by prior epilepsy
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("no prior epilepsy", "prior epilepsy"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by structural etiology
summary(pSERG[pSERG$structuraletiology==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.50 47.74 45.00 460.00
sd(pSERG[pSERG$structuraletiology==1, ]$BZDTIME.0)
## [1] 81.15584
summary(pSERG[pSERG$structuraletiology==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 17.00 71.78 45.00 1440.00
sd(pSERG[pSERG$structuraletiology==0, ]$BZDTIME.0)
## [1] 182.1918
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$structuraletiology 0.09882 1.10387 0.13280 0.744 0.457
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$structuraletiology 1.104 0.9059 0.8509 1.432
##
## Concordance= 0.509 (se = 0.016 )
## Rsquare= 0.002 (max possible= 1 )
## Likelihood ratio test= 0.55 on 1 df, p=0.5
## Wald test = 0.55 on 1 df, p=0.5
## Score (logrank) test = 0.55 on 1 df, p=0.5
# Figure time to first BZD by structural etiology
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("no structural etiology", "structural etiology"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by hospital onset
summary(pSERG[pSERG$HOSPITALONSET=="no", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.00 20.00 69.13 55.00 1264.00
sd(pSERG[pSERG$HOSPITALONSET=="no", ]$BZDTIME.0)
## [1] 155.9042
summary(pSERG[pSERG$HOSPITALONSET=="yes", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 57.53 25.00 1440.00
sd(pSERG[pSERG$HOSPITALONSET=="yes", ]$BZDTIME.0)
## [1] 176.3347
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$HOSPITALONSETyes 0.3299 1.3908 0.1260 2.619 0.00882 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$HOSPITALONSETyes 1.391 0.719 1.087 1.78
##
## Concordance= 0.563 (se = 0.016 )
## Rsquare= 0.022 (max possible= 1 )
## Likelihood ratio test= 6.56 on 1 df, p=0.01
## Wald test = 6.86 on 1 df, p=0.009
## Score (logrank) test = 6.92 on 1 df, p=0.009
# Figure time to first BZD by hospital onset
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("out of the hospital", "in the hospital"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by type of status epilepticus
summary(pSERG[pSERG$TYPESTATUS=="intermittent", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 20.00 85.71 60.00 1440.00
sd(pSERG[pSERG$TYPESTATUS=="intermittent", ]$BZDTIME.0)
## [1] 193.1872
summary(pSERG[pSERG$TYPESTATUS=="continuous", ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 6.00 14.00 23.31 28.00 180.00
sd(pSERG[pSERG$TYPESTATUS=="continuous", ]$BZDTIME.0)
## [1] 29.20201
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$TYPESTATUSintermittent -0.4991 0.6071 0.1296 -3.852 0.000117
##
## pSERG$TYPESTATUSintermittent ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$TYPESTATUSintermittent 0.6071 1.647 0.4709 0.7826
##
## Concordance= 0.536 (se = 0.017 )
## Rsquare= 0.046 (max possible= 1 )
## Likelihood ratio test= 14.12 on 1 df, p=2e-04
## Wald test = 14.84 on 1 df, p=1e-04
## Score (logrank) test = 15.13 on 1 df, p=1e-04
# Figure time to first BZD by type of status epilepticus
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("continuous", "intermittent"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by day versus night
summary(pSERG[pSERG$day==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 16.00 54.95 41.75 1132.00
sd(pSERG[pSERG$day==1, ]$BZDTIME.0)
## [1] 126.6443
summary(pSERG[pSERG$day==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 20.00 80.14 50.00 1440.00
sd(pSERG[pSERG$day==0, ]$BZDTIME.0)
## [1] 201.1023
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$day))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$day)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$day 0.1230 1.1308 0.1182 1.041 0.298
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$day 1.131 0.8843 0.8971 1.426
##
## Concordance= 0.509 (se = 0.018 )
## Rsquare= 0.004 (max possible= 1 )
## Likelihood ratio test= 1.09 on 1 df, p=0.3
## Wald test = 1.08 on 1 df, p=0.3
## Score (logrank) test = 1.08 on 1 df, p=0.3
# Figure time to first BZD by day versus night
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$day), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("night", "day"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by weekday versus weekend/holiday
summary(pSERG[pSERG$holiday==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 6.00 20.00 73.67 50.00 1264.00
sd(pSERG[pSERG$holiday==0, ]$BZDTIME.0)
## [1] 163.9659
summary(pSERG[pSERG$holiday==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 10.50 47.14 31.25 1440.00
sd(pSERG[pSERG$holiday==1, ]$BZDTIME.0)
## [1] 157.8114
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$holiday)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$holiday1 0.2851 1.3299 0.1265 2.254 0.0242 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$holiday1 1.33 0.7519 1.038 1.704
##
## Concordance= 0.541 (se = 0.016 )
## Rsquare= 0.016 (max possible= 1 )
## Likelihood ratio test= 4.89 on 1 df, p=0.03
## Wald test = 5.08 on 1 df, p=0.02
## Score (logrank) test = 5.11 on 1 df, p=0.02
# Figure time to first BZD by weekday versus weekend/holiday
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("weekday", "weekend/holiday"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
# Time to first treatment by season in the academic year
summary(pSERG[pSERG$earlyacademicyear==1, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 15.00 55.06 35.75 1264.00
sd(pSERG[pSERG$earlyacademicyear==1, ]$BZDTIME.0)
## [1] 140.9172
summary(pSERG[pSERG$earlyacademicyear==0, ]$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 6.75 20.00 74.70 55.00 1440.00
sd(pSERG[pSERG$earlyacademicyear==0, ]$BZDTIME.0)
## [1] 178.8734
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$earlyacademicyear 0.1752 1.1915 0.1162 1.508 0.131
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$earlyacademicyear 1.192 0.8393 0.9489 1.496
##
## Concordance= 0.529 (se = 0.018 )
## Rsquare= 0.008 (max possible= 1 )
## Likelihood ratio test= 2.26 on 1 df, p=0.1
## Wald test = 2.28 on 1 df, p=0.1
## Score (logrank) test = 2.28 on 1 df, p=0.1
# Figure time to first BZD by season in the academic year
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear), fun = "event",
conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3,
cex.axis = 1.3, cex.lab = 1.5,
frame.plot=FALSE,
xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("January-June", "July-December"),
col=c("violetred3", "seagreen3"),
lty=1,
lwd=3,
horiz=FALSE,
bty='n',
cex=1.3)
Multivariable model time to first BZD
# Initial model with the candidate variables with a p-value less than 0.25 on multivariable analysis:
#hospital onset, type of SE, holiday, and season in the academic year
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$earlyacademicyear))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$earlyacademicyear)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.0667 1.0690 0.1172 0.569 0.56938
## pSERG$HOSPITALONSETyes 0.3552 1.4265 0.1283 2.769 0.00562
## pSERG$TYPESTATUSintermittent -0.5736 0.5635 0.1320 -4.345 1.4e-05
## pSERG$holiday1 0.3075 1.3600 0.1296 2.372 0.01769
## pSERG$earlyacademicyear 0.1300 1.1388 0.1173 1.108 0.26793
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$earlyacademicyear
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0690 0.9355 0.8495 1.3451
## pSERG$HOSPITALONSETyes 1.4265 0.7010 1.1094 1.8342
## pSERG$TYPESTATUSintermittent 0.5635 1.7746 0.4350 0.7299
## pSERG$holiday1 1.3600 0.7353 1.0549 1.7535
## pSERG$earlyacademicyear 1.1388 0.8781 0.9049 1.4332
##
## Concordance= 0.601 (se = 0.021 )
## Rsquare= 0.096 (max possible= 1 )
## Likelihood ratio test= 30.28 on 5 df, p=1e-05
## Wald test = 31.37 on 5 df, p=8e-06
## Score (logrank) test = 31.82 on 5 df, p=6e-06
# Model without season in the academic year because p-value >0.1 in the multivariable model and it is not
#a confounder because it does not change the coefficient of intellectual disability by more than 20%
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.07959 1.08285 0.11669 0.682 0.49516
## pSERG$HOSPITALONSETyes 0.35566 1.42712 0.12855 2.767 0.00566
## pSERG$TYPESTATUSintermittent -0.58428 0.55751 0.13198 -4.427 9.56e-06
## pSERG$holiday1 0.31743 1.37360 0.12954 2.450 0.01427
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0828 0.9235 0.8615 1.3611
## pSERG$HOSPITALONSETyes 1.4271 0.7007 1.1093 1.8360
## pSERG$TYPESTATUSintermittent 0.5575 1.7937 0.4304 0.7221
## pSERG$holiday1 1.3736 0.7280 1.0656 1.7706
##
## Concordance= 0.599 (se = 0.021 )
## Rsquare= 0.092 (max possible= 1 )
## Likelihood ratio test= 29.06 on 4 df, p=8e-06
## Wald test = 30.03 on 4 df, p=5e-06
## Score (logrank) test = 30.5 on 4 df, p=4e-06
# Add non-candidate variables to the multivariable model: age (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$ageyears))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$ageyears)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.070733 1.073295 0.119648 0.591 0.55440
## pSERG$HOSPITALONSETyes 0.352977 1.423298 0.128753 2.741 0.00612
## pSERG$TYPESTATUSintermittent -0.591457 0.553520 0.133788 -4.421 9.83e-06
## pSERG$holiday1 0.317058 1.373082 0.129529 2.448 0.01437
## pSERG$ageyears 0.003761 1.003768 0.011263 0.334 0.73842
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$ageyears
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0733 0.9317 0.8489 1.3569
## pSERG$HOSPITALONSETyes 1.4233 0.7026 1.1059 1.8319
## pSERG$TYPESTATUSintermittent 0.5535 1.8066 0.4258 0.7195
## pSERG$holiday1 1.3731 0.7283 1.0652 1.7699
## pSERG$ageyears 1.0038 0.9962 0.9819 1.0262
##
## Concordance= 0.597 (se = 0.021 )
## Rsquare= 0.093 (max possible= 1 )
## Likelihood ratio test= 29.17 on 5 df, p=2e-05
## Wald test = 30.09 on 5 df, p=1e-05
## Score (logrank) test = 30.57 on 5 df, p=1e-05
# Add non-candidate variables to the multivariable model: sex (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$SEX))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$SEX)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.07884 1.08203 0.11666 0.676 0.49915
## pSERG$HOSPITALONSETyes 0.36019 1.43360 0.12885 2.795 0.00518
## pSERG$TYPESTATUSintermittent -0.57804 0.56100 0.13253 -4.362 1.29e-05
## pSERG$holiday1 0.32484 1.38380 0.13038 2.492 0.01272
## pSERG$SEXmale 0.05958 1.06139 0.12035 0.495 0.62054
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$SEXmale
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.082 0.9242 0.8609 1.3600
## pSERG$HOSPITALONSETyes 1.434 0.6975 1.1137 1.8455
## pSERG$TYPESTATUSintermittent 0.561 1.7825 0.4327 0.7274
## pSERG$holiday1 1.384 0.7226 1.0718 1.7867
## pSERG$SEXmale 1.061 0.9422 0.8384 1.3437
##
## Concordance= 0.599 (se = 0.021 )
## Rsquare= 0.093 (max possible= 1 )
## Likelihood ratio test= 29.3 on 5 df, p=2e-05
## Wald test = 30.3 on 5 df, p=1e-05
## Score (logrank) test = 30.76 on 5 df, p=1e-05
# Add non-candidate variables to the multivariable model: white (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$white))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$white)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.08515 1.08888 0.11695 0.728 0.46653
## pSERG$HOSPITALONSETyes 0.36022 1.43364 0.12877 2.797 0.00515
## pSERG$TYPESTATUSintermittent -0.57366 0.56346 0.13272 -4.322 1.54e-05
## pSERG$holiday1 0.31720 1.37328 0.12960 2.448 0.01438
## pSERG$white -0.08857 0.91524 0.12149 -0.729 0.46598
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$white
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0889 0.9184 0.8658 1.3694
## pSERG$HOSPITALONSETyes 1.4336 0.6975 1.1139 1.8452
## pSERG$TYPESTATUSintermittent 0.5635 1.7748 0.4344 0.7309
## pSERG$holiday1 1.3733 0.7282 1.0652 1.7704
## pSERG$white 0.9152 1.0926 0.7213 1.1613
##
## Concordance= 0.601 (se = 0.021 )
## Rsquare= 0.094 (max possible= 1 )
## Likelihood ratio test= 29.58 on 5 df, p=2e-05
## Wald test = 30.68 on 5 df, p=1e-05
## Score (logrank) test = 31.13 on 5 df, p=9e-06
# Add non-candidate variables to the multivariable model: prior epilepsy (enters the model: p-value>0.1 but confounder because it changes the coefficient for ID from 0.07959 to 0.03839, more than 20%)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$priorepilepsy)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.03839 1.03914 0.14013 0.274 0.78411
## pSERG$HOSPITALONSETyes 0.35980 1.43304 0.12878 2.794 0.00521
## pSERG$TYPESTATUSintermittent -0.59333 0.55249 0.13311 -4.457 8.3e-06
## pSERG$holiday1 0.31362 1.36837 0.12964 2.419 0.01556
## pSERG$priorepilepsy 0.07492 1.07780 0.14115 0.531 0.59554
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$priorepilepsy
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0391 0.9623 0.7896 1.3676
## pSERG$HOSPITALONSETyes 1.4330 0.6978 1.1134 1.8445
## pSERG$TYPESTATUSintermittent 0.5525 1.8100 0.4256 0.7172
## pSERG$holiday1 1.3684 0.7308 1.0613 1.7642
## pSERG$priorepilepsy 1.0778 0.9278 0.8173 1.4213
##
## Concordance= 0.598 (se = 0.021 )
## Rsquare= 0.093 (max possible= 1 )
## Likelihood ratio test= 29.34 on 5 df, p=2e-05
## Wald test = 30.22 on 5 df, p=1e-05
## Score (logrank) test = 30.72 on 5 df, p=1e-05
# Add non-candidate variables to the multivariable model: structural etiology (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$structuraletiology))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$priorepilepsy + pSERG$structuraletiology)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.04230 1.04320 0.13998 0.302 0.76253
## pSERG$HOSPITALONSETyes 0.34645 1.41404 0.12987 2.668 0.00764
## pSERG$TYPESTATUSintermittent -0.59226 0.55308 0.13285 -4.458 8.27e-06
## pSERG$holiday1 0.33052 1.39170 0.13144 2.515 0.01192
## pSERG$priorepilepsy 0.07784 1.08095 0.14087 0.553 0.58059
## pSERG$structuraletiology 0.11021 1.11652 0.13568 0.812 0.41660
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$priorepilepsy
## pSERG$structuraletiology
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0432 0.9586 0.7929 1.3725
## pSERG$HOSPITALONSETyes 1.4140 0.7072 1.0963 1.8239
## pSERG$TYPESTATUSintermittent 0.5531 1.8081 0.4263 0.7176
## pSERG$holiday1 1.3917 0.7185 1.0756 1.8006
## pSERG$priorepilepsy 1.0809 0.9251 0.8201 1.4247
## pSERG$structuraletiology 1.1165 0.8956 0.8558 1.4566
##
## Concordance= 0.598 (se = 0.021 )
## Rsquare= 0.095 (max possible= 1 )
## Likelihood ratio test= 29.99 on 6 df, p=4e-05
## Wald test = 30.97 on 6 df, p=3e-05
## Score (logrank) test = 31.42 on 6 df, p=2e-05
# Add non-candidate variables to the multivariable model: day versus night (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$day))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$priorepilepsy + pSERG$day)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.03354 1.03411 0.14046 0.239 0.81129
## pSERG$HOSPITALONSETyes 0.36118 1.43502 0.12846 2.812 0.00493
## pSERG$TYPESTATUSintermittent -0.58839 0.55522 0.13319 -4.418 9.98e-06
## pSERG$holiday1 0.31087 1.36461 0.12954 2.400 0.01641
## pSERG$priorepilepsy 0.09301 1.09747 0.14262 0.652 0.51431
## pSERG$day 0.11949 1.12692 0.11960 0.999 0.31774
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$priorepilepsy
## pSERG$day
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0341 0.9670 0.7852 1.3618
## pSERG$HOSPITALONSETyes 1.4350 0.6969 1.1156 1.8459
## pSERG$TYPESTATUSintermittent 0.5552 1.8011 0.4277 0.7208
## pSERG$holiday1 1.3646 0.7328 1.0586 1.7590
## pSERG$priorepilepsy 1.0975 0.9112 0.8298 1.4514
## pSERG$day 1.1269 0.8874 0.8914 1.4246
##
## Concordance= 0.6 (se = 0.021 )
## Rsquare= 0.096 (max possible= 1 )
## Likelihood ratio test= 30.34 on 6 df, p=3e-05
## Wald test = 31.4 on 6 df, p=2e-05
## Score (logrank) test = 31.88 on 6 df, p=2e-05
# Final multivariable model
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$priorepilepsy)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.03839 1.03914 0.14013 0.274 0.78411
## pSERG$HOSPITALONSETyes 0.35980 1.43304 0.12878 2.794 0.00521
## pSERG$TYPESTATUSintermittent -0.59333 0.55249 0.13311 -4.457 8.3e-06
## pSERG$holiday1 0.31362 1.36837 0.12964 2.419 0.01556
## pSERG$priorepilepsy 0.07492 1.07780 0.14115 0.531 0.59554
##
## pSERG$delay_or_cerebralpalsy
## pSERG$HOSPITALONSETyes **
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1 *
## pSERG$priorepilepsy
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0391 0.9623 0.7896 1.3676
## pSERG$HOSPITALONSETyes 1.4330 0.6978 1.1134 1.8445
## pSERG$TYPESTATUSintermittent 0.5525 1.8100 0.4256 0.7172
## pSERG$holiday1 1.3684 0.7308 1.0613 1.7642
## pSERG$priorepilepsy 1.0778 0.9278 0.8173 1.4213
##
## Concordance= 0.598 (se = 0.021 )
## Rsquare= 0.093 (max possible= 1 )
## Likelihood ratio test= 29.34 on 5 df, p=2e-05
## Wald test = 30.22 on 5 df, p=1e-05
## Score (logrank) test = 30.72 on 5 df, p=1e-05
# Test for proportionality
time.dep.zphfirstBZD <- cox.zph(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
time.dep.zphfirstBZD
## rho chisq p
## pSERG$delay_or_cerebralpalsy 0.0399 0.507 0.47628
## pSERG$HOSPITALONSETyes -0.1786 9.636 0.00191
## pSERG$TYPESTATUSintermittent -0.1269 5.220 0.02233
## pSERG$holiday1 -0.0249 0.189 0.66355
## pSERG$priorepilepsy -0.0433 0.609 0.43513
## GLOBAL NA 19.607 0.00148
# Proportional hazards assumtion for ID
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy), fun = "cloglog",
conf.int=FALSE, col = c("red", "blue"), lwd = 3,
cex.axis = 1.5, cex.lab = 1.5,
xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")
# Plot residual plots for ID
plot(time.dep.zphfirstBZD[1],
lwd = 2,
col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)
# Proportional hazards assumtion for hospital onset
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET), fun = "cloglog",
conf.int=FALSE, col = c("red", "blue"), lwd = 3,
cex.axis = 1.5, cex.lab = 1.5,
xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")
# Plot residual plots for hospital onset
plot(time.dep.zphfirstBZD[2],
lwd = 2,
col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)
# Proportional hazards assumtion for type of SE
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS), fun = "cloglog",
conf.int=FALSE, col = c("red", "blue"), lwd = 3,
cex.axis = 1.5, cex.lab = 1.5,
xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")
# Plot residual plots for type of SE
plot(time.dep.zphfirstBZD[3],
lwd = 2,
col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)
# Proportional hazards assumtion for holiday
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday), fun = "cloglog",
conf.int=FALSE, col = c("red", "blue"), lwd = 3,
cex.axis = 1.5, cex.lab = 1.5,
xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")
# Plot residual plots for holiday
plot(time.dep.zphfirstBZD[4],
lwd = 2,
col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)
# Proportional hazards assumtion for prior epilepsy
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy), fun = "cloglog",
conf.int=FALSE, col = c("red", "blue"), lwd = 3,
cex.axis = 1.5, cex.lab = 1.5,
xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")
# Plot residual plots for holiday
plot(time.dep.zphfirstBZD[5],
lwd = 2,
col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)
Multivariable model time to first BZD with interaction with hospital onset
# Create the variable HOSPITALONSETNUMERIC
pSERG$HOSPITALONSETNUMERIC <- ifelse(pSERG$HOSPITALONSET=="yes", 1, 0)
# Create the variable interactionIDhospitalonset
pSERG$interactionIDhospitalonset <- pSERG$delay_or_cerebralpalsy * pSERG$HOSPITALONSETNUMERIC
# Final multivariable model with interaction with hospital onset (non-significant)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$interactionIDhospitalonset))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy +
## pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday +
## pSERG$priorepilepsy + pSERG$interactionIDhospitalonset)
##
## n= 300, number of events= 300
##
## coef exp(coef) se(coef) z
## pSERG$delay_or_cerebralpalsy 0.01144 1.01150 0.15663 0.073
## pSERG$HOSPITALONSETyes 0.31039 1.36396 0.18259 1.700
## pSERG$TYPESTATUSintermittent -0.58623 0.55642 0.13441 -4.362
## pSERG$holiday1 0.31114 1.36498 0.12982 2.397
## pSERG$priorepilepsy 0.06861 1.07101 0.14205 0.483
## pSERG$interactionIDhospitalonset 0.09982 1.10497 0.25969 0.384
## Pr(>|z|)
## pSERG$delay_or_cerebralpalsy 0.9418
## pSERG$HOSPITALONSETyes 0.0891 .
## pSERG$TYPESTATUSintermittent 1.29e-05 ***
## pSERG$holiday1 0.0165 *
## pSERG$priorepilepsy 0.6291
## pSERG$interactionIDhospitalonset 0.7007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy 1.0115 0.9886 0.7441 1.3750
## pSERG$HOSPITALONSETyes 1.3640 0.7332 0.9536 1.9508
## pSERG$TYPESTATUSintermittent 0.5564 1.7972 0.4276 0.7241
## pSERG$holiday1 1.3650 0.7326 1.0583 1.7605
## pSERG$priorepilepsy 1.0710 0.9337 0.8107 1.4149
## pSERG$interactionIDhospitalonset 1.1050 0.9050 0.6642 1.8382
##
## Concordance= 0.596 (se = 0.021 )
## Rsquare= 0.094 (max possible= 1 )
## Likelihood ratio test= 29.49 on 6 df, p=5e-05
## Wald test = 30.51 on 6 df, p=3e-05
## Score (logrank) test = 31.07 on 6 df, p=2e-05